透過實作田納西-伊士曼製程(Tennessee Eastman Process, TEP)資料集,我們會帶大家了解一個機器學習專案是如何從零開始實作的,培養定義問題和解決問題的能力。在這個案例中,我們會學到:
1. 資料科學流程(Data Science Process)¶
2. 定義問題並了解資料來源¶
3. 撰寫程式前置作業¶
4. 資料讀取¶
5. 資料的初步觀察和前處理¶
6. 將資料建構成可以給模型的資料格式¶
7. 資料切分為訓練集和驗證集¶
8. 訓練模型:從簡單的模型開始¶
9. 評估剛剛建立好的模型¶
10. 如何改進結果¶
- 資料前處理-Feature scaling
- 模型訓練:選用 Tree-based 模型
- 模型訓練:選用 Ensemble 模型
- 資料特徵工程-Principle Component Analysis(PCA)
11. 異常檢測的特殊作法¶
- PCA 故障檢測
- Ensembling
12. 簡單的小總結¶
### Tennessee Eastman Process Simulation Dataset
### The diagram of process 整個化學反應製成主要有五個階段:反應器、冷凝器、循環壓縮機、分離器、汽提塔。主要反映過程如下:
| Variable number | Process variable | Type |
|---|---|---|
| IDV(1) | A/C feed ratio, B composition constant(stream 4) | Step |
| IDV(2) | B composition, A/C ratio constant(stream 4) | Step |
| IDV(3) | D feed temperature(stream 2) | Step |
| IDV(4) | Reactor cooling water inlet temperature | Step |
| IDV(5) | Condenser cooling water inlet temperature | Step |
| IDV(6) | A feed loss(stream 1) | Step |
| IDV(7) | C header pressure loss-reduced availability(stream 4) | Step |
| IDV(8) | A, B, C feed composition(stream 4) | Random variation |
| IDV(9) | D feed temperature(stream 2) | Random variation |
| IDV(10) | C feed temperature(stream 4) | Random variation |
| IDV(11) | Reactor cooling water inlet temperature | Random variation |
| IDV(12) | Condenser cooling water inlet temperature | Random variation |
| IDV(13) | Reaction kinetics | Slow drift |
| IDV(14) | Reactor cooling water valve | Sticking |
| IDV(15) | Condensor cooling water valve | Sticking |
| IDV(16) | Unknown | Unknown |
| IDV(17) | Unknown | Unknown |
| IDV(18) | Unknown | Unknown |
| IDV(19) | Unknown | Unknown |
| IDV(20) | Unknown | Unknown |
【程式用法】- import package
import numpy as np # 矩陣操作
import pandas as pd # 處理表格類型資料
import matplotlib.pyplot as plt # 視覺化
import seaborn as sns # 進階視覺化
np.set_printoptions(edgeitems=25, linewidth=150, formatter=dict(float=lambda x: "%.3g" % x)) # 設定顯示範圍
# TEP_FaultFree_training_100run.csv
!gdown --id 1ckTubYilJSW9q9cmp89NGpaarqQlt2NF
# TEP_Fault_training_25run.csv
!gdown --id 1ktdevLBIeFUSRwparVpyAEE7cTlKZxDg
# TEP_Fault_training_10run.csv
!gdown --id 1zwiqv6GsnCn3jbrHXZjRJhd9iW16Eeig
# TEP_Fault_testing_10run.csv
!gdown --id 1URWWeZQyO3FhS0U_nPuzMkqEb0EfX03Y
資料集為模擬運行,採樣時間均為每 3 分鐘一筆,故每小時將有 20 筆監測資料。訓練集每回運行時間為 25 小時;測試集則為 48 小時。其中蒐集異常資料時,訓練集會以正常操作開始值至第 1 小時模擬錯誤操作,因此前 20 筆為正常資料,而後 480 筆為異常資料;測試集會以正常操作開始值至第 8 小時模擬錯誤操作,因此前 160 筆為正常資料,而後 800 筆為異常資料
檔案說明:
【程式用法】- arguments of function
def function(X, a=1, b=None):
print(f'X={X}, a={a}, b={b}')
train_normal_df = pd.read_csv('TEP_FaultFree_training_100run.csv') # 將指定名稱的 csv 檔案讀入
train_fault_df = pd.read_csv('TEP_Fault_training_10run.csv')
train_normal_df.head(n=5) # 資料前 n 筆
train_normal_df.info() # 資料資訊:資料格式、筆數、欄位名稱及其形態等
train_fault_df.head(n=5)
xmeas_1~xmeas_41: 製程監控參數
| Description | Variable number | Base case value | Units |
|---|---|---|---|
| A feed(stream 1) | xmeas_1 | 0.25052 | kscmh |
| D feed(stream 2) | xmeas_2 | 3664.0 | kgh⁻¹ |
| E feed(stream 3) | xmeas_3 | 4509.3 | kgh⁻¹ |
| A and C feed(stream 4) | xmeas_4 | 9.3477 | kscmh |
| Recycle flow(stream 8) | xmeas_5 | 26.902 | kscmh |
| Reactor feed rate(stream 6) | xmeas_6 | 42.339 | kscmh |
| Reactor pressure | xmeas_7 | 2705.0 | kPa gauge |
| Reactor level | xmeas_8 | 75 | % |
| Reactor temperature | xmeas_9 | 120.4 | °C |
| Purge rate(stream 9) | xmeas_10 | 0.33712 | kscmh |
| Product separator temperature | xmeas_11 | 80.109 | °C |
| Product separator level | xmeas_12 | 50 | % |
| Product separator pressure | xmeas_13 | 2633.7 | kPa gauge |
| Product separator underflow(stream 10) | xmeas_14 | 25.16 | m³h⁻¹ |
| Stripper level | xmeas_15 | 50 | % |
| Stripper pressure | xmeas_16 | 3102.2 | kPa gauge |
| Stripper underflow(stream 11) | xmeas_17 | 22.949 | m³h⁻¹ |
| Stripper temperature | xmeas_18 | 65.731 | °C |
| Stripper steam flow | xmeas_19 | 230.31 | kgh⁻¹ |
| Compressor work | xmeas_20 | 341.43 | kW |
| Reactor cooling water outlet temperature | xmeas_21 | 94.599 | °C |
| Separator cooling water outlet temperature | xmeas_22 | 77.297 | °C |
| Reactor feed analysis(stream 6) | ||||
|---|---|---|---|---|
| Component | Variable number | Base case value | Units | Sampling frequency=0.1h |
| A | xmeas_23 | 32.188 | mol% | |
| B | xmeas_24 | 8.8933 | mol% | |
| C | xmeas_25 | 26.383 | mol% | |
| D | xmeas_26 | 26.383 | mol% | |
| E | xmeas_27 | 26.383 | mol% | |
| F | xmeas_28 | 26.383 | mol% | |
| Purge gas analysis(stream 9) | ||||
| Component | Variable number | Base case value | Units | Sampling frequency=0.1h |
| A | xmeas_29 | 32.958 | mol% | |
| B | xmeas_30 | 13.823 | mol% | |
| C | xmeas_31 | 23.978 | mol% | |
| D | xmeas_32 | 1.2565 | mol% | |
| E | xmeas_33 | 18.579 | mol% | |
| F | xmeas_34 | 2.2633 | mol% | |
| G | xmeas_35 | 4.8436 | mol% | |
| H | xmeas_36 | 2.2986 | mol% | |
| Product analysis(stream 11) | ||||
| Component | Variable number | Base case value | Units | Sampling frequency=0.25h |
| D | xmeas_37 | 0.01787 | mol% | |
| E | xmeas_38 | 0.83570 | mol% | |
| F | xmeas_39 | 0.09858 | mol% | |
| G | xmeas_40 | 53.724 | mol% | |
| H | xmeas_41 | 43.828 | mol% |
| Description | Variable number | Base case value(%) | Low limit | High limit | Units |
|---|---|---|---|---|---|
| D feed flow(stream 2) | xmv_1 | 63.053 | 0 | 5811 | kgh⁻¹ |
| E feed flow(stream 3) | xmv_2 | 53.980 | 0 | 8354 | kgh⁻¹ |
| A feed flow(stream 1) | xmv_3 | 24.644 | 0 | 1.017 | kscmh |
| A and C feed flow(stream 4) | xmv_4 | 61.302 | 0 | 15.25 | kscmh |
| Compressor recycle valve | xmv_5 | 22.210 | 0 | 100 | % |
| Purge valve(stream 9) | xmv_6 | 40.064 | 0 | 100 | % |
| Separator pot liquid flow(stream 10) | xmv_7 | 38.100 | 0 | 65.71 | m³h⁻¹ |
| Stripper liquid product flow(stream 11) | xmv_8 | 46.534 | 0 | 49.10 | m³h⁻¹ |
| Stripper steam valve | xmv_9 | 47.446 | 0 | 100 | % |
| Reactor cooling water flow | xmv_10 | 41.106 | 0 | 227.1 | m³h⁻¹ |
| Condenser cooling water flow | xmv_11 | 18.114 | 0 | 272.6 | m³h⁻¹ |
對於此資料集 X 以及 y 分別是:
train_normal_df.describe() # 資料描述:預設針對數值型態欄位統計相關數據,包含計數、平均值、標準差、最小、第一四分位距、第二四分位距、第三四分位距、最大值等。
TEP_FaultFree_training_100run.csv 僅存在 faultNumber=0(Normal),而各個監測值的值域範圍差異極大。另外,不存在著欄位中僅有一種值。
train_fault_df.describe() # 資料描述:預設針對數值型態欄位統計相關數據,包含計數、平均值、標準差、最小、第一四分位距、第二四分位距、第三四分位距、最大值等。
TEP_Fault_training_10run.csv 中 faultNumber=1~20(Process Faults)
取一小部分的資料來做比較:Normal data v.s. Fault data
【程式用法】- slicing
a = list(range(0, 15))
print('a:\t ', a)
print('a[:]\t ', a[:])
print('a[:2]:\t ', a[:2])
print('a[2:6]:\t ', a[2:6])
print('a[1:]:\t ', a[1:])
print('a[2:8:2]:', a[2:8:2])
【程式用法】- pd.DataFrame
【程式用法】 - conditional statement
輸出為 boolean (布林值):True / False
結合 array 索引 → 只取用為 True 的索引資料
a = np.arange(10)
print('a:\t', a)
print('a>5:\t', a>5)
print('a[a>5]:\t', a[a>5])
sample_train_normal = train_normal_df[train_normal_df.simulationRun==1] # 透過條件限制 simulationRun=1 取用部分資料
sample_train_fault = train_fault_df[train_fault_df.simulationRun==1]
【程式用法】- for loop
for idx, i in enumerate(range(6, 12)):
print(f'{idx+1} th loop: {i}')
plt.figure(figsize=(15, 20)) # 設定畫布大小
for idx, i in enumerate(range(0, 5)): # 查看第 i 個特徵的分布
plt.subplot(5, 1, idx+1) # 在 5 列 1 行的畫布上,選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(sample_train_normal['sample'],
sample_train_normal.iloc[:, i+3],
label='Normal') # 繪製折線圖
for j in range(2): # 查看第 j 種 process fault
plt.plot(sample_train_fault.loc[sample_train_fault.faultNumber==j+1, 'sample'],
sample_train_fault.loc[sample_train_fault.faultNumber==j+1, f'xmeas_{i+1}'],
label=f'Fault_{j+1}') # 繪製折線圖
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.title(sample_train_normal.columns[i+3]) # 設定標題
plt.tight_layout() # 自動保持子圖之間的間距
plt.figure(figsize=(15, 20)) # 設定畫布大小
for idx, i in enumerate(range(0, 5)): # 查看第 i 個特徵的分布
plt.subplot(5, 1, idx+1) # 在 5 列 1 行的畫布上,選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(sample_train_normal['sample'],
sample_train_normal.iloc[:, i+3],
label='Normal') # 繪製折線圖
for j in range(2, 4):
plt.plot(sample_train_fault.loc[sample_train_fault.faultNumber==(j+1), 'sample'],
sample_train_fault.loc[sample_train_fault.faultNumber==(j+1), f'xmeas_{i+1}'],
label=f'Fault_{j+1}') # 繪製折線圖
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.title(sample_train_normal.columns[i+3]) # 設定標題
plt.tight_layout() # 自動保持子圖之間的間距
先觀察少部分的資料:對於 faultNumber = 3 and 4 與 Normal 的前 5 個特徵分布十分相近,可以再進一步觀察其他特徵的分布差異並分析此分類是否存在問題。
comb_df = pd.concat((sample_train_normal.iloc[20:, :], sample_train_fault[sample_train_fault['sample']>=20])) # 合併 normal 和 fault data
comb_df['faultNumber'] = comb_df['faultNumber'].astype('int') # 將欄位 faultNumber 轉換為 int(整數)型態
plt.figure(figsize=(20, 12)) # 設定畫布大小
for idx, i in enumerate(range(3)): # 查看第 i 個特徵分布
plt.subplot(3, 1, idx+1) # 在 3 列 1 行的畫布上,選擇編號為 idx+1 的位置
sns.boxplot(data=comb_df[['faultNumber', f'xmeas_{i+1}']],
x='faultNumber',
y=f'xmeas_{i+1}') # 繪製盒鬚圖
plt.tight_layout() # 自動保持子圖之間的間距
plt.figure(figsize=(25, 12)) # 設定畫布大小
for idx, i in enumerate(range(3)): # 查看第 i 個特徵分布
plt.subplot(3, 1, idx+1) # 在 3 列 1 行的畫布上,選擇編號為 idx+1 的位置
sns.violinplot(data=comb_df[['faultNumber', f'xmeas_{i+1}']],
x='faultNumber',
y=f'xmeas_{i+1}') # 繪製小提琴圖
plt.tight_layout() # 自動保持子圖之間的間距
corr = comb_df.iloc[:,[0]+[i for i in range(3, 55)]].corr(method='kendall') # 計算 faultNumber 以及 xmeas, xmv 間的相關係數
plt.figure(figsize=(30, 20)) # 設定畫布大小
sns.heatmap(corr) # 繪製熱圖(heatmap)
comb_df['faultNumber'] = comb_df['faultNumber'].astype('str') # 將欄位 faultNumber 轉換為 str(字串)型態
new_comb_df = pd.get_dummies(comb_df, columns=['faultNumber']) # 將 faultNumber 拆成 20 個欄位(One-hot encoding:屬於該類別,給值 1 反之為 0)
new_comb_df.head()
corr = new_comb_df.iloc[:,2:].corr(method='kendall') # 計算各種 faultNumber 以及 xmeas, xmv 間的相關係數
plt.figure(figsize=(30, 20)) # 繪製畫布大小
sns.heatmap(corr) # 繪製熱圖
train_all_10run = pd.concat((train_normal_df[train_normal_df['simulationRun']<=10], train_fault_df[train_fault_df['sample']>20])) # 合併 normal 和 fault data
train_all_10run.describe()
from sklearn.model_selection import train_test_split
X_train, X_valid, y_train, y_valid = train_test_split(train_all_10run.iloc[:, 3:].values, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17) # 隨機方式
#, stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
print(f'The shape of X_train: {X_train.shape}\t y_train: {y_train.shape}')
print(f'The shape of X_valid: {X_train.shape}\t y_valid: {y_valid.shape}')
unique, counts = np.unique(train_all_10run['faultNumber'].values, return_counts=True) # 計算各種 fault 有多少筆
plt.bar(unique, counts) # 繪製直條圖
unique, counts = np.unique(y_train, return_counts=True) # 計算各種 fault 有多少筆
plt.bar(unique, counts) # 繪製直條圖
【程式用法】- Sklearn model
from sklearn.linear_model import LogisticRegression
model = LogisticRegression() # 建立模型
model.fit(X_train, y_train) # 訓練模型
from sklearn.metrics import mean_squared_error, classification_report, confusion_matrix
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
# model.score(X_valid, y_valid)
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣(Confusion matrix)
先看訓練集,再看驗證集!
from sklearn.preprocessing import StandardScaler # (X-mean(X))/std(X)
from sklearn.preprocessing import MinMaxScaler # (X-min(X))/(max(X)-min(X))
def data_preprocessing(df_input, train=True, sc=None):
if train: # 若是 training data,則按照輸入資料做 feature scaling
sc = StandardScaler()
# sc = MinMaxScaler()
df = sc.fit_transform(df_input)
else: # 若是 testing data,則按照 training scale 做 feature scaling
df = sc.transform(df_input)
return df, sc
train_all_10run.describe()
train_transform_data, sc = data_preprocessing(train_all_10run.iloc[:, 3:]) # feature scaling for 52 features
pd.DataFrame(train_transform_data, columns=train_all_10run.columns[3:]).describe()
X_train, X_valid, y_train, y_valid = train_test_split(train_transform_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
model = LogisticRegression() # 建立模型
model.fit(X_train, y_train) # 訓練模型
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣(Confusion matrix)
from sklearn.tree import DecisionTreeClassifier
X_train, X_valid, y_train, y_valid = train_test_split(train_transform_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
model = DecisionTreeClassifier() # 建立模型
model.fit(X_train, y_train) # 訓練模型
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣(Confusion matrix)
模型過度擬合在訓練資料集上(綠色線)達到很高的準確率,導致在未參與訓練的驗證資料集上表現不佳
抑制 Overfitting 的方法:
from lightgbm import LGBMClassifier
X_train, X_valid, y_train, y_valid = train_test_split(train_transform_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
model = LGBMClassifier() # 建立模型
model.fit(X_train, y_train) # 訓練模型
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣(Confusion matrix)
from sklearn.decomposition import PCA
pca = PCA(n_components=0.99) # n_components: int, float or ‘mle’
pca.fit(train_transform_data) # 計算共變異數矩陣並進行特徵分解
pca_transformed_data = pca.transform(train_transform_data) # pca 降維轉換
pca_transformed_data.shape
pca.explained_variance_ratio_
pca.explained_variance_ratio_.sum()
X_train, X_valid, y_train, y_valid = train_test_split(pca_transformed_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
model = LGBMClassifier() # 建立模型
model.fit(X_train, y_train) # 訓練模型
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣(Confusion matrix)
通常在做異常檢測時,蒐集到的資料集為大量的正常資料,以及少數的異常資料。因此,若以單純的分類做法,往往會因為資料不平均而導致未能檢測出異常資料。
然而在異常檢測問題上,會利用以下的想法,來訓練模型:
為了更貼近實際狀況-正常資料與異常資料比例差異懸殊,將會使用模擬 100 回正常訓練資料,以及 10 回異常資料來做演示
train_normal_df.describe()
sc = StandardScaler()
sc_normal = sc.fit_transform(train_normal_df.iloc[:, 3:]) # 按照 normal data scale 做 feature scaling
pca_normal = PCA(n_components=0.99) # n_components: int, float or ‘mle’
transformed_normal_data = pca_normal.fit_transform(sc_normal) # 計算共變異數矩陣並進行特徵分解後,做降維轉換
pca_normal.explained_variance_ratio_
transformed_normal_data.shape
T2_normal_data = (transformed_normal_data**2 / pca_normal.explained_variance_).sum(axis=1) # 計算 T^2 統計量
sc_fault = sc.transform(train_fault_df.iloc[:, 3:]) # 按照 normal data scale 對 fault data 做 feature scaling
transformed_fault_data = pca_normal.transform(sc_fault) # 按照 normal data 的降維方式對 fault data 做降維
transformed_fault_data.shape
T2_fault_data = (transformed_fault_data**2 / pca_normal.explained_variance_).sum(axis=1) # 計算 T^2 統計量
plt.figure(figsize=(25, 40)) # 設定畫布大小
for idx, i in enumerate(range(20)): # 查看第 i 種 fault 的 T^2 分布
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上,選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(range(len(T2_normal_data[:500])), # 繪製折線圖
T2_normal_data[:500],
label='Normal')
plt.plot(range(len(T2_fault_data[i*500:(i+1)*500])), # 繪製折線圖
T2_fault_data[i*500:(i+1)*500],
label=f'Fault_{i+1}')
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
normal_project = np.einsum('ik,kf->ifk',
transformed_normal_data,
pca_normal.components_).sum(axis=2) # 計算 normal data 的投影向量
Q_normal_data = ((normal_project-(train_normal_df.iloc[:, 3:]-train_normal_df.mean()[3:]))**2).sum(axis=1) # 計算 Q 統計量
fault_project = np.einsum('ik,kf->ifk',
transformed_fault_data,
pca_normal.components_).sum(axis=2) # 計算 fault data 的投影向量
Q_fault_data = ((fault_project-(train_fault_df.iloc[:, 3:]-train_fault_df.mean()[3:]))**2).sum(axis=1) # 計算 Q 統計量
plt.figure(figsize=(25, 40)) # 設定畫布大小
for idx, i in enumerate(range(20)): # 查看第 i 種 fault 的 Q 分布
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上,選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(range(len(Q_normal_data[:500])), # 繪製折線圖
Q_normal_data[:500],
label='Normal')
plt.plot(range(len(Q_fault_data[i*500:(i+1)*500])), # 繪製折線圖
Q_fault_data[i*500:(i+1)*500],
label=f'Fault_{i+1}')
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
sc = StandardScaler()
sc_normal_df = pd.DataFrame(sc.fit_transform(train_normal_df.iloc[:, 3:]), # 按照 normal data scale 做 feature scaling
columns=train_normal_df.iloc[:, 3:].columns)
import tqdm
from lightgbm import LGBMRegressor
def train(df, cols_to_predict):
models = {} # 各個模型存放位置
pbar = tqdm.tqdm_notebook(cols_to_predict)
for col in pbar:
pbar.set_description(f'Training model for {col}')
model = LGBMRegressor(learning_rate=0.1) # 建立模型
tr_x = df.drop([col],axis=1) # 取其中 N-1 個特徵作為 X
target = df[col] # 取剩餘的 1 個特徵作為 y
model.fit(X=tr_x, y=target) # 訓練模型
models[col] = model # 將訓練好的模型存放至 models
return models
def predict(models, df, cols_to_predict):
preds = []
for col in cols_to_predict:
test_x = df.drop([col], axis=1) # 取其中 N-1 個特徵作為 X
pred = models[col].predict(test_x) # 預測剩餘的 1 個特徵
preds.append(pred)
return preds
sc_normal_df.describe()
features_to_predict = sc_normal_df.columns
models = train(sc_normal_df, features_to_predict) # 透過自定義的 train 訓練模型
def get_mse(sample, preds):
return np.square((sample.loc[:,features_to_predict] - np.transpose(preds))**2).mean(axis=1) # 計算 mean square error
plt.figure(figsize=(25,40)) # 設定畫布大小
plt.title('MSE for a normal and fault samples') # 設定標題
normal_Run = np.random.randint(100)+1 # 隨機取正常資料的其中一回合
print(f'Normal simulationRun = {normal_Run}')
normal_sample = pd.DataFrame(sc.transform(train_normal_df[train_normal_df.simulationRun==normal_Run].iloc[:, 3:]),
columns=train_normal_df.columns[3:])
normal_preds = predict(models, normal_sample, features_to_predict) # 透過自定義的 predict 預測結果
fault_Run = np.random.randint(10)+1 # 隨機取異常資料的其中一回合
print(f'fault simulationRun = {fault_Run}')
for idx, i in enumerate(range(20)):
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上,選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(get_mse(normal_sample, normal_preds),
label='Normal') # 得到預測結果與原始特徵之間的 MSE 值
faulty_sample = pd.DataFrame(sc.transform(train_fault_df[(train_fault_df.simulationRun==fault_Run) & (train_fault_df.faultNumber==(i+1))].iloc[:,3:]),
columns=train_fault_df.columns[3:])
faulty_preds = predict(models, faulty_sample, features_to_predict) # 透過自定義的 predict 預測結果
plt.plot(get_mse(faulty_sample, faulty_preds),
label=f'Fault_{i+1}') # 得到預測結果與原始特徵之間的 MSE 值
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
plt.figure(figsize=(25,40)) # 設定畫布大小
plt.title('MSE for a normal and fault samples') # 設定標題
normal_Run = np.random.randint(100)+1 # 隨機取正常資料的其中一回合
print(f'Normal simulationRun = {normal_Run}')
normal_sample = pd.DataFrame(sc.transform(train_normal_df[train_normal_df.simulationRun==normal_Run].iloc[:, 3:]),
columns=train_normal_df.columns[3:])
normal_preds = predict(models, normal_sample, features_to_predict) # 透過自定義的 predict 預測結果
fault_Run = np.random.randint(10)+1 # 隨機取異常資料的其中一回合
print(f'fault simulationRun = {fault_Run}')
for idx, i in enumerate(range(20)):
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上,選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(get_mse(normal_sample, normal_preds),
label='Normal') # 得到預測結果與原始特徵之間的 MSE 值
faulty_sample = pd.DataFrame(sc.transform(train_fault_df[(train_fault_df.simulationRun==fault_Run) & (train_fault_df.faultNumber==(i+1))].iloc[:,3:]),
columns=train_fault_df.columns[3:])
faulty_preds = predict(models, faulty_sample, features_to_predict) # 透過自定義的 predict 預測結果
plt.plot(get_mse(faulty_sample, faulty_preds),
label=f'Fault_{i+1}') # 得到預測結果與原始特徵之間的 MSE 值
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
此資料集也關乎時間上的關聯性,由 Fault 17 所呈現的圖形明顯看出具有週期性的變化,因此在機器學習的做法當中,可以加入往前推算一段時間的所有特徵,或者時間相關的因素在特徵當中,例如:過往時間段的平均值、最大值、最小值等等。
然而在深度學習的做法當中,有專門為時間序列所設計的 RNN 模型,也是未來可以模型方面嘗試的方向。
簡而言之,決定模型性能的常常是資料的好壞。而資料的好壞又會被資料處理的方式影響。能否找到好的資料處理方式既切合商業目標,並讓模型好學習是資料科學專案的重點。